library(sf)
library(mapview)
library(dplyr)Importing additional spatial data
Mined Units
The mined units can be read in directly with st_read.
mined <- st_read('https://ca.dep.state.fl.us/arcgis/rest/services/OpenData/MMP_MINEDUNITS/MapServer/0/query?outFields=*&where=1%3D1&f=geojson')Reading layer `OGRGeoJSON' from data source
`https://ca.dep.state.fl.us/arcgis/rest/services/OpenData/MMP_MINEDUNITS/MapServer/0/query?outFields=*&where=1%3D1&f=geojson'
using driver `GeoJSON'
Simple feature collection with 371 features and 10 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -82.95369 ymin: 27.46063 xmax: -81.69582 ymax: 30.51925
Geodetic CRS: WGS 84mapview(mined)Brownfield sites
Download the zipped kml file to a temporary directory.
# url with zipped kml
urlin <- 'https://ordsext.epa.gov/FLA/www3/acres_frs.kmz'
# download file
tmp1 <- tempfile(fileext = ".kmz")
download.file(url = urlin, destfile = tmp1, method = 'curl')Unzip the kmz file.
tmp2 <- tempdir()
unzip(tmp1, exdir = tmp2)Get the name of the kml file to read.
lyr <- unzip(tmp1, list = T)$Name
fl <- paste(c(tmp2, lyr), collapse = "\\")
fl <- gsub('\\\\', '/', fl)Read the kml file with st_read and drop Z dimension with st_zm. Here, the Tampa sites are loaded. You can view all possible locations in the kml file with st_layers.
dat <- st_read(fl, layer = 'TAMPA') %>%
st_zm()Reading layer `TAMPA' from data source `/tmp/RtmpqdPSoL/ACRES_FRS.KML' using driver `LIBKML'
Simple feature collection with 103 features and 12 fields
Geometry type: POINT
Dimension: XYZ
Bounding box: xmin: -82.51857 ymin: 27.87108 xmax: -82.3533 ymax: 28.07826
z_range: zmin: 0 zmax: 0
Geodetic CRS: WGS 84mapview(dat)To import all layers in the kml file, identify the layer names and loop through them to add to a single object. This takes several hours. The data are saved as an .Rdata object and .csv file for later use. The data include only the site name and location in decimal degrees.
# layer names
alllyr <- st_layers(fl)$name
strt <- Sys.time()
out <- NULL
for(i in alllyr){
# counter
cat(i, which(i == alllyr), 'of', length(alllyr), '\n')
print(Sys.time() - strt)
# import each layer
dat <- st_read(fl, i, quiet = T)[, c('Name')]
# append to same object
out <- rbind(out, dat)
}
# save as RData object
allbfld <- out %>% st_zm()
save(allbfld, file = 'data/allbfld.RData', compress = 'xz')
# save as csv
allbfldcsv <- allbfld %>%
mutate(
lon = st_coordinates(.)[,1],
lat = st_coordinates(.)[,2]
) %>%
st_set_geometry(NULL)
write.csv(allbfldcsv, 'data/allbfldcsv.csv', row.names = F)View all brownfield sites.
load(file = 'data/allbfld.RData')
mapview(allbfld, legend = F, col.regions = 'grey')